APPROACH TO A STATIONARY STATE IN AN 

EXTERNAL FIELD 



A.Alastuey* J.Piasecki** 
* ENS Lyon, CNRS, France 
**Institute of Theoretical Physics, 
University of Warsaw, Hoza 69, 
00 681 Warsaw, Poland 
(Dated: March 4, 2010) 

Abstract 

We study relaxation towards a stationary out of equilibrium state by analyzing a one- 
dimensional stochastic process followed by a particle accelerated by an external field and 
propagating through a thermal bath. The effect of collisions is described within Boltz- 
mann's kinetic theory. We present analytical solutions for the Maxwell gas and for the 
very hard particle model. The exponentially fast relaxation of the velocity distribution 
toward the stationary form is demonstrated. In the reference frame moving with constant 
drift velocity the hydrodynamic diffusive mode is shown to govern the distribution in the 
position space. We show that the exact value of the diffusion coefficient for any value of 
the field is correctly predicted by the Green-Kubo autocorrelation formula generalized to 
the stationary state. 

PACS numbers: 05.30.-d, 05.70.Ce, 52.25.Kn 
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I. INTRODUCTION 



The present paper is devoted to the the study of a stochastic process followed by a 
particle moving through a scattering thermal bath while accelerated by an external 
field. The field prevents the particle from acquiring the Maxwell distribution of 
the bath. Our aim here is not only to establish the precise form of the stationary 
velocity distribution, as it was e.g. the case in the analysis presented in jlj, but also 
to answer the physically relevant question of the dynamics of approach towards the 
long-time asymptotic state. The evolution of the distribution in position space will 
be thus also discussed. 

We consider a one-dimensional dynamics described by the Boltzmann kinetic 
equation 

/ d d d \ f 

\di + V dr + a dv ) ^ V] ^ = ^ P I dw \ v ~ f( r > w; ^ ^ ~ v; ^ ^ w "> ^ 

(1) 

Here f(r, v; t) is the probability density for finding the propagating particle at point r 
with velocity v at time t. The thermal bath particles are not coupled to the external 
field. Before binary encounters with the accelerated particle they are assumed to be 
in an equilibrium state with uniform temperature T and density p 
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Here 4>{v) is the Maxwell distribution, and 
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denotes the corresponding thermal velocity. The differential operator on the left- 
hand side of ([T]) generates motion with a constant acceleration a. The accelerated 
motion is permanently perturbed by instantaneous exchanges of velocities with ther- 
malized bath particles. This is modeled by the Boltzmann collision term on the right 
hand side of equation ([1]), which accounts for elastic encounters between equal mass 
particles. The collision frequency depends therein on the absolute relative velocity 
\v — w\ through a simple power law with exponent 7. Finally v int is some character- 
istic velocity of the underlying interparticle interaction. 



In the case of hard rods (7 = 1) the factor \v — w\ is the main source of difficulties 
in the attemps to rigorously determine the evolution of f(r, v; t), since it prevents the 
effective use of Laplace and Fourier transformations. It was thus quite remarkable 
that a stationary velocity distribution could be analytically determined in that case, 
leading in particu>ar to an explicit expression for the current at any valne of the 

Jj. In that case, kinetic equation (CQ) has been solved exactly 



external acceleration 



only at zero temperature where 4>(v)\t=o = 5(v) Also, when 0(t>) is replaced 
by the distribution [S(v — t>o) + S(v + Vq)]/2 with a discrete velocity spectrum ±t> , 
an explicit analytic solution has been derived and analyzed in 3f and J4]. The 
physically relevant conclusions from those works can be summarized as follows 

(i) the approach to the asymptotic stationary velocity distribution is exponen- 
tially fast 

(ii) in the reference system moving with average velocity, the hydrodynamic dif- 
fusion mode governs the spreading of the distribution in position space 

(iii) the Green-Kubo autocorrelation formula for the diffusion coefficient applies in 
the non-equilibrium steady state 

Our aim is to show that the general features (i)-(iii) persist when <ft(v) is the 
Maxwell distribution with temperature T > 0. However, in the present study, we 
restrict the analysis to cases 7 = and 7 = 2, which are much simpler than the 
hard-rod one. Indeed, it turns out that the Fourier-Laplace transformation can 
then be effectively used to solve the initial value problem for equation ([1]). The 
simplifications occuring when 7 = or 7 = 2 have been already exploited in other 
studies: for recent applications to granular fluids, see e.g. [sJ-Ql and references 
quoted therein. 

In terms of dimensionless variables 

w = v/v th , x = rp(v th /v int ) 7_1 , r = tpv th (v th /v int y~ 1 , (4) 

the kinetic equation (TjQ) takes the form 
/ d d d \ f 

(— — \~ w ~q — \~ e ~Q— J F( x i w 'j t ) = / du\w — m| 7 [F(x, u; r)$(w) — F(x, w; t)$(m)] , 

(5) 



where $(«;) is the dimensionless normalized gaussian 



Hw) = ^=e~ w2 l\ (6) 



and e is the dimensionless parameter 

e = {v th v int ) - (7) 

proportional to the ratio between the energy amp' 1 provided to the particle on a 
mean free path, and thermal energy ksT. That parameter can thus be looked upon 
as a measure of the strength of the field. Integration of ([5]) over the position space 
yields the kinetic equation for the velocity distribution 



G(w;r) = J dxF(x, w; r) , 

which reads 

The paper is organized as follows. In Section II, we consider the so-called Maxwell 
gas (7 = 0). The explicit solution of the kinetic equation (J5J) enables a thorough 
discussion of the approach to the stationary state, together with a study of the 
structure of the stationary velocity distribution. In Section III, we proceed to a 
similar analysis for the very hard particle model (7 = 2). Section IV contains 
conclusions. Some calculations have been relegated to Appendices. 

II. THE MAXWELL GAS 

We consider here the simple version 7 = of equation (jSj). One usually then refers 
to the Maxwell gas dynamics, in which the collision frequency does not depend on 
the speed of approach (see e.g. pi]). This case can be viewed upon as a very crude 
approximation to the hard rod dynamics (7 = 1) obtained by replacing the relative 
speed \v — c\ of colliding particles by constant thermal velocity v th , while v int is 
identified with v th . Here, kinetic equation (jSJ) takes the form 
d d d \ f 

dr + W dx + € &w J W]T ^ = j & U [ F ( X > M ? t)$(w) - F(x, w; r)$(tt)] 

= M (x;r)$(w) - F(x,w;t) (9) 



where M (x; r) denotes the zeroth moment 

M (x;t)= / duF(x, u; r) . (10) 



Equation ([9]) can be conveniently rewritten as an integral equation 

F(x, w; t) = e~ T F(x — wt + er 2 /2, w — er; 0) 

+ / dr]e~ r '§(w - er])M (x - wrj + er] 2 /2; r - rf) , (11) 
Jo 

with an explicit dependence on the initial condition F(x, w; 0). Integration of equa- 
tion ( TTTT) over x yields 

G(w; r) = J dxF(x, w; r) = e~ T G in (w - er) + N J dr/e^$(w - er]) , (12) 

where G in (w) = G(w;Q) is the initial condition, and N = J dw J dxF(x,w,T) = 
J dwG(w; t) is the conserved normalization factor. 



A. Stationary solution and relaxation of the velocity distribution 

Putting iVo = 1 in formula ffl2l yields the evolution law for the normalized 
velocity distribution 

G(w;t)= [ dxF(x, w; r) = e~ T G in (w - er) + / dr]e' v ^(w - erf) , (13) 



The first term on the right hand side of ( fl3l) describes the decaying memory of the 
initial distribution : G in (w) propagates in the direction of the field with constant 
velocity e, while its amplitude is exponentially damped. Clearly, for times r 3> 1 
that term can be neglected. 

The second term in formula ([TBI describes the approach to the asymptotic sta- 
tionary distribution 



G st (w) = G(w; oo) 



2e 

where 
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G at {w) ~ - exp ( — - - ) (15) 



FIG. 1: Stationary velocity distribution G st (w) for three values of e. 

is the familiar error function. It is interesting to compare the decay-law of G st (w) 
at large velocities, to that corresponding to the case of hard-rod collisions. Using 
expression (|T4|) we find the asymptotic formula 

1 f I w s 

7 exp U*~7 

when w — > +oo. In contradistinction to the hard-rod case governed by an e- 
dependent gaussian law (see l|) we find here a purely exponential decay. The 
thermal bath is unable to impose via collisions its own gaussian decay because of 
insufficient collision frequency. The replacement of the relative speed in the Boltz- 
mann collision operator by thermal velocity implies thus qualitative changes in the 
shape of the stationary velocity distribution. The plot of G st (w) for different values 
of e is shown in Fig. [TJ 

Basic properties (i)-(iii) discussed in the Introduction turn out to be valid. In- 
deed, the inequality 

G st (w) - / di] e~ v §(w - er]) = / d •qe~" n $(u> - er]) < (16) 

Jo Jt e 

displays an uniform exponentially fast approach towards the stationary state. In 
particular, using formula (|T3|) . we find that the average velocity < w > (r) ap- 
proaches the asymptotic value 

< w > Bt = e (17) 
6 



according to 

J dw w G(w; t) = e + e~ r [< w > in — e] 



< W > T 



We encounter here an exceptional situation where the linear response is exact for 
any value of the external field. 

Equation (1T21 with N put equal to zero can be used for the evaluation of the 
time-displaced velocity autocorrelation function 

r(r) =< [w(f)- < w >JH0)- < w > st ] > st . (19) 

where < ... > st denotes the average over stationary state ( fl4l) . The calculation 
presented in Appendix |B] provides the formula 

r(r) = e" T [l + e 2 ] , (20) 

which yields a remarkably simple field dependence of the diffusion coefficient 

POO 

D(e)= / drr(r) = l + e 2 . (21) 



Jo 

B. Relaxation of density: appearence of a hydrodynamic mode 

Let us turn now to the analysis of the evolution of the normalized density 
n(x; t) = M (x; r) in position space. It turns out that one can solve the complete in- 
tegral equation ( II ip by applying to both sides Fourier and Laplace transformations. 
If we set 

poo I* 

F(k,w;z) = J dTe~ ZT J dxe~ lkx F(x,w;t), (22) 

we find 

f°° r ( t 2 \ 

F(k,w;z) = J dr exp — ik ( wr — e— J — (z + l)r 

{F in {k,w-eT)+n{k;z)${w-eT)} , (23) 

where n(k; z) is the Fourier-Laplace transform of n(x; r), and 
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denotes the spatial Fourier transform of the initial condition. Equation ( 12 3 j) when 
integrated over the velocity space yields the formula 

i r . r° . r / r 2 \ 1 * 

-zA; I wr — e— 1 — (z + 1)t io — er) 

(24) 
with 



n(/c; 2) = ^jT^T j dw I dr exp 







/■oo 

z) = 1 — dr exp 



(z + l)r - (ike + £; 2 )y 



(25) 



The insertion of ( )24|) into ( 123|) provides a complete solution for F(k,w;z) corre- 
sponding to a given initial condition. 



Formula (1241) shows that the time-dependence of the spatial distribution is defined 
by roots of the function ((k; z). In order to find the long-time hydrodynamic mode 
Zhy(fc)) we have to look for the root of ((k; z) which approaches when k — Y 0. If 
we assume the asymptotic form 

^hy(^) = C\k + c 2 k 2 + o{k 2 ) when k — > , 

we find a unique self-consistent solution to equation ((k; z) = of the form 

Zhy (k) = -iek - (1 + e 2 )k 2 + o{k 2 ) = -iek - D{e)k 2 + o{k 2 ) . (26) 

It has the structure of a propagating diffusive mode. It is important to note that the 
diffusion coefficient D(e) equals (1 + e 2 ) in accordance with the Green-Kubo result 
(|2~TT) . We thus see that, in the reference system moving with constant velocity e, a 
classical diffusion process takes place in position space. 

It has been argued in the literature that, in general, z\ iy (k) is not an analytic 
function of k at k = (see e.g. Ref. [9(). Here, that question can be precisely 
investigated as follows. According to the integral expression ( 1251) of £(k;z), the 
hydrodynamic mode is a function of £ = ike + k 2 . By combining differentiations 
with respect to £ under the integral sign with integration by parts, we find that 
z hy(0 satisfies the second order differential equation 

,^£hy = l dz^ 
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Then, since Zh y (0) = 0, we find that ^hy(C) can be formally represented by an infinite 
entire series in £, 

oo 

MO = E c ^ n ' (28) 

n=l 

with c\ = — 1, c 2 = 1 and 

|c„+i| > 2™" 1 n! for n > 2 . 

Thus, the radius of convergence of Taylor series (12"5|) is zero, so £ = is a singular 
point of function Zh y (0> as we ^ as = is a singular point of function z^j(k). The 
nature of that singularity can be found by rewriting the root equation defining Zh y (£) 
as the implicit equation 

The introduction of function a/£ requires to define cut-lines ending at points k = 
and k = —ie which are the two roots of equation = 0. Since the integral in 
the r.h.s. of expression ( 125|) diverges for k imaginary of the form k = iq with q > 
or q < — e, it is natural to define such cut-lines as [z'0,zoo[ and ] — zoo, —ie]. The 
corresponding choice of determination for \/£ is defined by a/£(& + ) = i\Jqt + q 2 for 
k + = + + iq with q > 0, where \Jqe + q 2 is the usual real positive square root of 
the real positive number (qe + q 2 ). Notice that, when complex variable k makes a 
complete tour around point k = starting from k + = + + iq on one side of the 
cut-line and ending at k~ = CT + iq on the other side (with vanishing difference 
k + — k~), v£pO changes sign from to — — with obvious notations. 
As shown by adding both implicit equations fl29|) for A; + and k~ respectively, z^ y 
does not reduce to z^ y . The difference (z£ y — z^ y ) is of order exp(— l/(2|fc|e)), so 
k = is an essential singularity. 

III. VERY HARD PARTICLES 

Another interesting case is that of the so-called very hard particle model, where 
the collision frequency is proportional to the kinetic energy of the relative motion 
of the colliding pair. The corresponding exponent in the collision term of the Boltz- 
mann equation ([[]) is now 7 = 2. This allows us to simplify the resolution of the 

9 



kinetic equation. Owing to this fact, the very hard particle model, similarly to the 
Maxwell gas, has been studied in numerous works (see e.g. [l(J-[llj], and references 
given therein). 

Using dimensionless variables @, we thus write the kinetic equation as 

F(x,w;r) — J du\w — u\ 2 [F(x, u; r)$(w) — F(x, w; r)$(?i)] 

(30) 

= [it? M (x; t) - 2wM l (x; t) + M 2 {x; t)]$(w) - {w 2 + l)F(x, w; T) 
where the moments Mj(x; r) (J = 1, 2, ...) are defined by 

M j (x;r)= dww j F(x,w;r) . (31) 



d_ +w _^ + e _d_\ 
dr dx dw J 



The evolution equation of the velocity distribution G(w; r) becomes 

— + e— ) G(w; r) = [N 2 {r) - 2wN 1 (r) + w 2 N ]$(w) - (w 2 + l)G{w; r) , (32) 



s dr dw t 
with the integrated moments 



Nj(r) = Jdx M 3 [x- r), j = 0, 1, 2 . (33) 

Notice that the integrated zeroth moment does not depend on time since the evolu- 
tion conserves the initial normalization condition 



N (t) = j dw j dxF(x,w,T) = N . 

Hence, when F(x,w; r) is a normalized probability density N Q (r) = N = 1. 

The simplification related to the choice 7 = 2, and more generally when 7 is an 
even integer, concerns the collision term in the general kinetic equation (pQ) which 
can be expressed in such cases in terms of a finite number of moments of the dis- 
tribution function. The resolution of that equation becomes then straightforward 
within standard methods (see Appendix lAf) . 

A. Laplace transform of the velocity distribution 

The expression for the Laplace transform of the normalized velocity distribution 
follows directly from the general formula (1A4I) derived in Appendix [A] by putting 
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k = 0, and choosing M (0,z) = Nq(z) = 1/z. Within definition 



S(w; z) = (z + l)w + 



(34) 



for the function S(k,w; z) evaluated at k = (see definition flA2j) ). we find 



eG{w; z) = -®(w) + / du exp{[S{u; z) - S{w; z))/e) {G in {u) 

Z J -oo 

+ [N 2 (z) - 2uN 1 (z) + {€U ~ Z z ~ l) ]<$>{u)} . 
The two functions N\{z) and N 2 (z) satisfy the system of equations 



(35) 



= A${0; z) + [N 2 (z) - (z + l)/z]A 00 (0; z) + [e/z - 2iV 1 (z)]A 01 (0; z) 
eNx(z) = ^'(O; z) + [N 2 (z) - (z + l)/z]A 10 (0; z) + [e/z - 2N 1 (z)]A n (0\ ^36) 

which is identical to (1A8I) taken at k — 0, while 

/f*W 
dw du exp{[S(u;z) - S(w; z)}/e} w j u l $(u) . (37) 
J — oo 

Analogous formula holds for (0; z) with the Maxwell distribution <3>(u) replaced 
by the initial condition G(u; 0) = G in (u). Once system ( |36|) has been solved, the 
insertion of the resulting expressions for Nx(z) and N 2 (z) into formula (1351) yields 
eventually an explicit solution of the kinetic equation for the velocity distribution 



G(w; z) = -A_> + - 
z e 



du exp {[^(it; z) — S{w\ z)]/e} 

i 

x {GM + [A t (z) u - B e (z))$(u)} . (38) 



With the shorthand notations Aji(z) = Aji(0;z) and Aj^'(z) = Aji(Q;z), the for- 
mulae for coefficients A e (z) and B e (z) read 

r 2 



AJz) 



A(z) 



A 00 (z) - 2A 00 (z)A^(z) + 2A 10 (z)A^(z) 



and 



BAz) 



A(z) 



A 01 (z) + eA^(z) + 2A u (z)A^(z) - 2A 01 (z)A^(z 



(39) 



(40) 



where A(z), in accordance with the definition given in HAlOj) . is 



A{z) = eA 00 (z) + 2 {A 00 {z)A u {z) - A 10 {z)A 01 {z)) . 



(41) 
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B. Stationary solution 



At large times, r — > oo, we expect the velocity distribution to reach some sta- 
tionary state G st (w) = G(w;oo). This can be easily checked by investigating the 
behaviour of G(w; z) in the neighbourhood of z = at fixed velocity w. 

All integrals over u in formula (|38|) do converge for any complex value of z. 
Moreover, all their derivatives with respect to z are also well defined, as shown by 
differentiation under the integral sign. Thus, such integrals are entire functions of 
z. The sole quantities in expression fl38|) which become singular at z = are the 
coefficients A e (z) and B e (z), and obviously the term <$>(w)/z. In fact, both A e (z) 
and B e (z) exhibit simple poles at z = 0. Hence, the stationary solution of the kinetic 
equation (1321) does emerge when r — > oo, and it is given by the residue of the simple 
pole of G(w; z) at z — 0, namely 

~S(u;0) - S(w;0)' 



G st (w) = $H + — J du exp 



[A 00 (0)u -Aoi(0)Mu). 

(42) 

In that expression, Ajj(0) and A(0) are the non-zero values at z = of the analytic 
functions Aij(z) = Ay(0; z) and A(z) = A(0; z). Formula rtlSjl does not depend on 
initial condition G in . All initial conditions evolve towards the same unique stationary 
distribution (142]) . It can be checked that the direct resolution of the static version 
of kinetic equation fl32|) obtained by setting dG/dr = does provide formula fT42|) . 

Since the external field accelerates the particle, the stationary solution is asym- 
metric with respect to the reflection w — > —w, and positive velocities are favoured. 
This leads to a finite current 

(w) st = J dwwG st {w) = -^-r[A 00 {0)A 11 (0) - Aoi(0)Aio(0)] . (43) 

The asymptotic expansion at large velocities of G st (w), inferred from formula f|42|) . 
reads 



V Z7T 



1 + !^) +0( i) 



when \w\ — > oo . (44) 



A(0)w v w 2 ' 

Therefore, the external field does not influence the leading large- velocity behaviour 
of G st (w), which is identical to that of the thermal bath. Its effects only arise in the 

12 




FIG. 2: Stationary velocity distribution G at (w) for three values of e. 

first correction to the leading behaviour which is smaller by a factor of order 1/w. 
The stationary distribution is drawn in Fig. [2] for several increasing field strengths, 
e = 1, e = 10 and e = 100. 



Let us study now the limit e — > which corresponds to a weak external field. 
The main contributions to the integrals over u in (I38I) arise from the region close to 
w. That observation motivates the use of a new integration variable y = (w — u)/e. 
The Taylor expansions of the resulting integrands in powers of e generate then entire 
series in e, the first terms of which read 



du u $(«) exp 



S{u; 0) -S{w; 0) 



and 



du exp 



S(u,0) - S(w,0) 



e— ^ +0(e 

1 + w z 



1 + U7 2 



+ 0(6^) 



(45) 



(46) 



Consequently, also functions Aij(0) and A(0) can be represented by power series in 
e as they are obtained by calculating appropriate moments of expansions ( 14"5]) and 
( |46|) over the velocity space. The corresponding small-e expansion of the stationary 
velocity distribution reads 

b w 



GJw) =$(w)+e 



1 + w 2 



(47) 
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where 



w 2 



1 + 2 / dw- -${w) 

1 + w z 



Of course, at e = 0, G Bt (w) reduces to the Maxwell distribution. The first correction 
is of order e, as expected from linear response theory. The corresponding current 
Q4"2l) reduces to 

(w) st = ae + 0(e 2 ), (48) 
where the conductivity a is given by 

a =i(l -b) (49) 

It will be shown in the sequel that a = Dq = D(e = 0), where D(e) is the diffusion 
coefficient given by the Green-Kubo formula. 

Consider now the strong field limit e — > oo. The corresponding behaviours 
of Aij(0) and A(0) are derived from the integral representations obtained in Ap- 
pendix We then find at fixed w 



£ -l/3 rw 

G st (w) = jcc- — — / dw$(w)exp 



S{u,0) - S(w,0) 



+ 0(e- 2 / 3 ) (50) 



For w of order 1, the dominant term in the large-e expansion of the integral in ( 15 Op 
reduces to 

dri*(ti) = i(l + Erf^ 

and thus varies from to 1 around the origin w — 0. For larger values of the velocity, 
w ~ e 1 / 3 , that integral behaves as exp(— w 3 /(3e). The next term in the expansion 
(J5QD remains of order e -2 / 3 . Thus, when e — > oo at fixed e~ 1//3 u> the stationary 
solution is given by 

GM - 0(w) roo - 6 ^ 3 exp [-(6-^)3/3] , (51) 
J a V exp^ y I o) 

where 9 is the Heaviside step function. The whole distribution is shifted toward 
high velocities w ~ e 1 / 3 , so that the resulting current (|43p is of the same order of 
magnitude, i.e. 

(«>>- ~ 3l r ( i ( /3) 3) eV3 when e ^ 00 ' ( 52 ) 

where T is the Euler Gamma function. That behavior can be recovered within 
the following simple interpretation. At strong fields, the average velocity of the 
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FIG. 3: Average current (w) B t as a function of e. The dashed line represents the linear 
Kubo term in the small-e expansion (|48p with conductivity a ~ 0.2039. The dotted line 
describes asymptotic formula §2$ with 3 1 / 3 r(2/3)/r(l/3) ~ 0.7290 valid in the limit 



particle becomes large compared to the thermal velocity of scatterers. Since at 
each collision the particle exchanges its velocity with a thermalized scatterer, the 
variation of particle velocity between two successive collisions is of the order of (v) st . 
On the other hand, in the stationary state the same velocity variation is due to the 
acceleration a coming from the external field, so it is of the order ar coll where r coll 
is the mean time between two successive collisions. This time can be reasonably 
estimated as the inverse collision frequency for a relative velocity \v — c| of order 
(v) st . The consistency of those estimations requires the relation 



which indeed implies the e 1//3 -behaviour (|5"2|) of the average velocity in dimensionless 
units. Contrary to the Maxwell case where the current remains linear in the applied 
field, here the current deviates from its linear-response form when the field increases 
: it grows more slowly because collisions are more efficient in dissipating the energy 
input of the field. In Fig. [3J, we plot (w) st as a function of e. 



e — > oo. 




(53) 
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C. Relaxation towards the stationary solution 



Let us study now the relaxation of the velocity distribution G(w; r) towards 
the stationary solution G st (w). The decay of [G(w,r) — G st (w)] when r — » oo is 
controlled by the singularities of G(w; z) in the complex plane, different from the 
pole at z = 0. As already mentioned, all integrals in expression (1381) are entire 
functions of z, so the singularities at z ^ arise only in the coefficients A t (z) and 
B e (z). Thus, the first important conclusion is that the relaxation is uniform for the 
whole velocity spectrum. 

According to expressions (1391) and ( HU1) defining A e (z) and B e (z) respectively, the 
singularities of those coefficients at points z 7^ 0, correspond to zeros of the function 
A (2) given by expression (jUJ). Since the analytic functions Aij(z) and A(z) do not 
depend on initial condition G in . the relaxation is an intrinsic dynamical process, as 
expected. 

After some algebra detailed in Appendix EJ we find that A(z) reduces to the 
Laplace transform 

POO 

A(z) = e 2 dyf £ (y)ex P (-zy) (54) 



of the real, positive, and monotonously decreasing function 

2,, v f 2 (l+3y) / 3 yH2 + v)\ 

f M) = (l + »)(l + 2»)^ 6XP ~ £ 6(TWJ ' <55) 

Owing to the fast decay of f e (y) the integral fl54l) converges for any z, so A (2) is 
an entire function of z. Also, the monotonic decay of f e (y) and its positivity imply 
some general properties for the roots of A(z). First of all, A(z) cannot vanish 
for > 0. Moreover, as A(z) is strictly positive for 2; real, the zeros of A(z) 

appear in complex conjugate pairs, while they are isolated with strictly negative real 
parts and nonvanishing imaginary parts. Consequently, the long-time relaxation of 
the velocity distribution is governed by the pair of zeros which is closest to the 
imaginary axis. Noting them as z ± = — A ± iu with oj ^ and < A, we conclude 
that G(w; t) relaxes towards G st {w) via exponentially damped oscillations 

G(w; t) — G 3t (w) ~ G(w) cos[o;r + v( w )] ex P( — ^ r ) ? when r — > 00 (56) 
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where C(w) and r)(w) are an amplitude and a phase respectively. It should be 
noticed that both functions C{w) and r](w) depend on initial conditions. 

At a given value of e, the zeros z ± are found by solving numerically the equation 
A{z ) =0. In the weak- or strong-field limits, we can derive asymptotic formulae 
for such zeros as follows. First, as indicated by numerically computing z^ for small 
values of e, z ± collapse to z = —1 when e — > 0. The corresponding asymptotical 
behaviour can be derived by noting that, for z close to zo, the leading contributions 
to A(z) in integral ( )54l) arise from large values of y. Then, we set y = £/e 2//3 and 
z = — 1 + s e 2 / 3 , which provide 

o ,5/3 poo 

A(-l + se 2 / 3 )~| d£r 1/2 exp(- S £-e7i2) (57) 

when e — > at fixed s. By numerical methods, we find the pair of complex conju- 
gated zeros Sq of integral 



POD 

/ d£r 1/2 exp(- S £-£7i2) 



which are the closest to the imaginary axis. Therefore, when e — > 0, damping factor 
A(e) goes to 1 according to 

A(e) = 1 - e 2 / 3 + o(e 2 / 3 ) (58) 

with 9£(sq) ~ —1.169, while frequency u(e) vanishes as S(s^) e 2//3 with S(s^) ~ 
2.026. Notice that for fixed z, not located on the real half-axis ] — oo, —1], A(z) 
behaves as 

A(z)~e 2 A (z) (59) 

when e — > 0, with 



" (*+l)/2 



X 



l-Erf( v /(^ + l)/2^ 
V / 2^(7Tl)e (2+1)/2 (l - Erf U(z + l)/2\\j . (60) 



Here, \/ (z + l)/2 is defined as the usual real positive square root y/(x~+ l)/2 for 
real z = x belonging to the half axis x > — 1, while the complementary half-axis 
z = x<— lisa cut-line ending at the branching point z = —1. That point is the 
singular point of 1/A (z) closest to the imaginary axis, as strongly suggested by a 
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numerical search of the zeros of A (z). Therefore, both A(e) and oj(e) are continuous 
functions of e at e = with A(0) = 1 and uj(0) = 0. At e = 0, the exponentially 
damped oscillating decay (J56"j) becomes an exponentially damped monotonic decay 
multiplied by power- law t~ 3 / 2 . That power-law arises from the presence of a singular 
term of order \/(z + l)/2 in the expansion of G(w; z) around the branching point 
z = -1. 

When e — > oo, the zeros of A(z) are obtained by simultaneously changing y to 
£/e 2 / 3 in the integral and by rescaling z as e 2 / 3 s. This provides 

A(e 2/3 s) ~ e 4/3 A oc (s) when e -> oo at fixed s, (61) 

with 

/■oo 

A 00 (s)=/ d£exp(- S £-£ 3 /3) . (62) 

Therefore, when e — )■ oo, z 1 * 1 behave as 2; ± ~ e 2//3 s^, where s^, are the zeros of A^s) 
closest to the imaginary axis. The corresponding large-e asymptotical behaviour of 
the damping factor A(e) is 

A(e) = -^)e 2 / 3 + o(e 2 / 3 ) (63) 

with 3?(s^) — —2.726, while frequency u(e) diverges as 3(s+ ) e 2 / 3 with 3(s+ ) ~ 
6.260. Notice that the relaxation time A _1 (e) goes to zero as e _2//3 , like the average 
time between collisions r coll ~ (v) st /a used in our simple heuristic derivation of the 
e-dependence of the stationary current in the strong field limit. In Fig. HI we draw 
the damping factor A(e) function of e. 

D. Relaxation of density in position space 

In Appendix |X] we derive an explicit formula for the zeroth moment M (k; z) of 
the distribution F(k,w;z) which contains all information on the evolution of the 
spatial density of the propagating particle. The formula flA9j) clearly reveals the 
presence of a hydrodynamic pole in M (k; z), namely the root of equation 

z + (k 2 + iek) U{k; z) = (64) 
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Small-E behaviour 



FIG. 4: Damping factor A(e) as a function of e. The dashed and dotted lines represent 
the asymptotical behaviours ([58]) and (j63|) at small and large e respectively. 



where 



U(k;z) 



A u (k; z)A 00 (k; z) - A 10 (k; z)A 01 (k; z) 

eA 00 (k; z) + 2[A n (k; z)A 00 (k; z) - A 10 (k; z)A 01 (k; z)\ 



(65) 



If we consider the small-/c limit and if we assume the asymptotic form 

z hy (k) = -ick - D(e) k 2 + 0(k 2 ) 



(66) 



for the hydrodynamic root, we find immediately from equation (I64p the formula 

c = el7(0;0). (67) 

This shows that the mode propagates with the average stationary velocity (w) st = 
eU(0] 0) derived in expression f|43|) . 



In order to infer the formula for the diffusion coefficient D(e), it is necessary 
to calculate the term linear in variable k in the expansion of function U (k; z) at 
z = —ick. Indeed, equation (164")) implies the equality 

D(e) = U(0; 0) + le ^rU(k; -ick)\ k=0 . (68) 

Taking into account the structure (165]) of U (k; z) we find the formula 

D( £ ) = ^ + A °° ][A ' llAm ~ A ' mAl ^ + A Mo^ ~ A '^ A °^ ( 6 9) 



A 2 
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where all Aji and A are taken at k = z = 0, and where 

A' jh = te^A jt (k; -ick)\ k=0 . (70) 

A particularly useful representation of the derivative appearing in expression (I70p 
can be deduced from formulae ( 1A2I) and ( 1A6I) defining functions Aji(k; z). An inte- 
gration by parts yields 

A'= dw du{u-c) dvw j v l exp{[S(0,v;0) - S(Q,w;0)]/e}$(v) . (71) 



It is quite remarkable that equation (171 j) allows us to establish a relation between 
the quantities A'- t and the stationary velocity distribution G st (w). Indeed, using 
equation (j42l . we readily obtain the equalities 



dw / duexp{[S(0,v;) - 5(0, w; 0)]/e}(u - c) G st (u) 



A i - cA 00 + ^[A 00 A' Q1 - A 01 A' Q0 ] = Joi (72) 



and 



/w 
du exp{[S'(0, v; ) — 5(0, w] 0)]/e} w (u — c) G Bt (u) 
-oo 



= A n - cA 10 + ^[Ax,^ - A 01 A' 10 ] = J n . (73) 

Then, we find that the linear combination (A 00 Ju — A w J i) of integrals Jn and Joi 
reduces to 

A n A 00 - A W A 01 + - {AoofAooA'n - A 01 A' W ] - A W [A 00 A' 01 - A 01 A' 00 ]} . (74) 



The comparison of that expression with equation ( 1691) leads to the compact final 
result 

n/ x ^00^11—^10^01 /„ c s 

£>( e ) = ^ • (75) 

The above formula involves, via coefficients Jn and J i, averages over the stationary 
velocity distribution. In fact, we show in Appendix [B] that expression ( 1751) follows 
by extending, to the present out-of-equilibrium stationary state, the familiar Green- 
Kubo relation between the diffusion coefficient and the velocity fluctuations. That 
important fact is one of the main observations of the present study. 
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When e — > 0, the behaviour of D(e) is easily infered by inserting the small-e 
expansion ( )47|) of the stationary velocity distribution G st (w) into formula ( 1751) . We 
find that -D(e) goes to conductivity cr f|49|) as quoted above, with a negative e 2 - 
correction. When e — )■ oo, we can use the large-e form (l5Tj) of G st (u>) for evaluating 
coefficients Jn and J i- Using also the corresponding behaviours of coefficients A 00 
and A w , we eventually obtain that -D(e) goes to the finite value 

= r3(i/3)-9r ( i/3)r(2/3) + 6r»(2/3) 

2r 3 (i/3) v ; 

The external field dependence of the diffusion coefficient -D(e) is shown in Fig. |5j 

The expansion ( l66l) of -2hy(^) can be pursued beyond the fc 2 -diffusion term, by 
expanding function U (k; z) in double entire series with respect to z and k. According 
to the integral expression of functions Aji(k; z) derived in AppendixEU all coefficients 
of those double series are finite. This implies that the hydrodynamic root z^ y {k) of 
equation ( 164)) can be formally represented by an entire series in k, namely 

CO 

Zhy(k) = ^a n k n , 

n=l 

with «i = — ic and «2 = — D(e). Coefficient a n (n > 3) can be straightforwardly 
computed once lowest-order coefficients a p with 1 < p < n — 1 have been determined. 
As shown by that calculation, all coefficients are obviously finite. Therefore, and 
similarly to what happens in the Maxwell case, only positive integer powers of k 
appear in the small-fc expansion of Zhy(&)- Now, we are not able to determine the 
radius of convergence of that expansion, so we cannot conclude about the analyticity 
of function zt y (k). However, we notice that, contrarily to the Maxwell case, the 
integrals defining Aji(k; z) remain well-defined for any complex value of k, as soon 
as e 7^ (see Appendic [C]). This suggests that z^y(k) might be an analytic function 
of k at k = 0, except for e = 0, in which case k = should be a singular point. 

IV. CONCLUDING COMMENTS 

The idea of this work was to perform a detailed study of the approach to an out- 
of-equilibrium stationary state, by considering systems for which analytic solutions 
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FIG. 5: Diffusion coefficient D{e) as a function of e. The dotted line represents the 
constant asymptotic value D^. 

can be derived. To this end we solved, within Boltzmann's kinetic theory, the 
one- dimensional initial value problem for the distribution of a particle accelerated 
by a constant external field and suffering elastic collisions with thermalized bath 
particles. Our exact results for the Maxwell model and for the very hard particle 
model support the general picture mentioned in the Introduction: 

• a uniform exponentially fast relaxation of the velocity distribution 

• diffusive spreading in space in the reference system moving with stationary 
flow 

• equality between the diffusion coefficient appearing in the hydrodynamic mode 
and the one given by the generalized Green-Kubo formula 

Although both models display the same phenomena listed above, the variations 
of the respective quantities of interest with respect to e are different. First we notice 
that, as far as deformations of the equilibrium Maxwell distribution are concerned, 
the external field is much less efficient for very hard particles. This is well illustrated 
by comparing figures [T] and [2] : for the Maxwell system, a significative deformation 
of $ is found for e = 5, while for the very-hard particle model a similar deformation 
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is observed for e = 100. This can be easily interpreted as follows. The collision 
frequency for very hard particles becomes much larger than its Maxwell gas coun- 
terpart when the external field increases, so it costs more energy to maintain a 
stationary distribution far from the equilibrium one. That mechanism also explains 
various related observations. For instance, the large- velocity behaviour of G at (w) 
is identical to the equilibrium Gaussian for very hard particles, while it takes an 
exponential form in the Maxwell gas. Also, the average current (w) st increases more 
slowly when e — > oo for very hard particles, and the corresponding relaxation time 
A _1 (e) vanishes instead of remaining constant for the Maxwell gas. 

Among the above phenomena, the emergence of a symmetric diffusion process in 
the moving reference frame is quite remarkable. In such a frame, there is some kind 
of cancellation between the action of the external field and the effects of collisions 
induced by the counterflow of bath particles with velocity u* ath = — (v) st . The 
corresponding diffusion coefficient D(e) increases whith e for the Maxwell gas (case 
7 = 0), while it decreases and saturates to a finite value for very hard particles (case 
7 = 2). Therefore, beyond the previous cancellation, it seems that the large number 
of collisions for 7 = 2 shrink equilibrium fluctuations. On the contrary, for 7 = 0, 
since D(e) diverges when e — > 00, the residual effect of collisions in the reference 
frame seems to vanish and particles tend to behave as if they were free. 

We expect that the same qualitative picture should be valid in the hard rod case 
which corresponds to the intermediate value 7 = 1 of the exponent 7 in equation (CQ). 
The quantitative behaviours should interpolate between those described for 7 = 
and 7 = 2. For instance, the stationary distribution G st (w) computed in Ref. []| 
displays a large- velocity asymptotic behaviour which is indeed intermediary between 
those derived here for 7 = and 7 = 2. Also, the average current < v > st is of order 
e 1 / 2 for e large, which lies between the e- and e 1//3 -behaviours found for 7 = and 
7 = 2 respectively. Notice that the e 1,/3 -behaviour for 7 = 2 can be retrieved within 
a selfconsistent argument, which uses in an essential way the existence of the velocity 
scale related to the particle-particle interaction. Whereas the thermal velocity scale 
becomes irrelevant when e — > 00, the interaction scale remains important. In the 
case of hard rods such an interaction scale does not show in the kinetic equation, 
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and the unique combination of parameters having the dimension of velocity is y/a/ p, 
which does provide a different strong field behaviour of < v > st with order e 1//2 . 

Appendix A: Solution of the kinetic equation for very hard particles 

Applying to equation (15Uj) Fourier and Laplace transformations, we find 

d ~ 

e——F(k, w; z) + (z + 1 + ikw + w 2 )F(k, w; z) 
aw 

= F in (k, w) + [M 2 (k; z) - 2wM 1 (k; z) + w 2 M (k; z))${w) (Al) 

where Mj(k; z) is the double Fourier-Laplace transform of the j f -moment Mj(x; r) 
defined in expression fl3Tl) . while F in (k, w) is the spatial Fourier transform of the ini- 
tial condition F in (x, w) = F(x, w; 0). The first order equation (lAlj) can be rewritten 
in an integral form with the use of function 

1 w 2 

S(k,w;z) =w(z + l) + -w 3 + ik— , (A2) 

namely 

/w 
du exp{[S'(fc, u; z) — S(k 1 w; z)]/e} {F in (k, u) 
-oo 

+ [M 2 (k; z) - 2uM 1 (k; z) + u 2 M (k; z)]$(u)} . (A3) 



Using then the relation 
d 



e— — h eu — (z + 1) — iku — u 2 
du 



exp{S(k,u;z)/e}$(u) = 



to evaluate the term involving u 2 in the right hand side of equation (1A3j) . we even- 
tually find the more convenient integral equation 



eF(k,w;z) = eM (h; z)$(w) + / du exp{[S(k, u; z) - S(k, w; z)]/e} 

J — oo 

x {F in (k, u) + [M 2 {k; z) - 2uM 1 (k; z) + (eu -iku-z- l)M {k; z)]$(u)} . (A4) 
Equation ( 1A4j) has to be considered together with the continuity equation 

zM (k;z) + ikM 1 (k;z) = M (k;0) . (A5) 
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In order to determine the unknown functions Mj(k;z) (j = 0,1,2), we complete 
equation (1A5j) with the zeroth and the first moments of equation (1A4I) . In the 
resulting system of linear equations, the integrals 

dw duw j u l exp{[S(k,u;z)-S(k,w;z)]/e}^(u) (A6) 

J — oo 

and 

dw duw J u l exp{[S(k,u;z) - S{k,w; z)]/e}F in {k,u) (A7) 

J — oo 

appear. That system reads 

= + [M 2 - (z + l)M }A 00 + [(e - ik)M - 2M 1 ]A 01 (A8) 
eM x = A ( (q + [M 2 - (z + 1)M ]A 10 + [(e - ik)M - 2Mi]An 

The explicit solution for the zeroth moment reads 

M = AM (k;0)-tk[A^A 00 -A 10 A^] 
zA + (ki + zek)(A n A 00 -A 10 A 01 ) 1 ' 

where 

A = eA 00 + 2(A U A 00 - A 10 ^oi) • (A10) 

The formula for the first moment follows directly from the continuity equation flA5j) . 
and then M 2 can be derived directly from (1A8j) . The insertion of the formulae for 
the first three moments into the relation flA4j) yields the complete solution for the 
distribution F(k, w; z) for any initial condition. 

Appendix B: Evaluation of the diffusion coefficient via Green-Kubo theory 

1. Velocity autocorrelation function of the Maxwell gas 

In order to evaluate the velocity autocorrelation function, we use the integral 
representation 

r(r) =< [w(t)- < w > st )[w(0)- < w > 8t ] > st = J dw w G(w; r) (Bl) 

where G(w; r) is the solution of kinetic equation f|T2|) corresponding to the initial 
condition 

G in (w) = (w- < w > Bt )G st (w) = (w- e)G st (w) (B2) 
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Here we find 

N = J dw G{w; t) = J dw G{w; 0) = , (B3) 
so that equation (fl2"|) takes a particularly simple form 

G(w; r) = e- T G in (w - er) = e~ T (w - e)G st (w - er) . (B4) 

Using then the explicit form (Till) of G st , we obtain 

/rOC 
dwwe' T (w-e)J drje^^w - erj) = e~ T (l + e 2 ) , (B5) 

which leads to the simple formula for the diffusion coefficient 

POO 

D(e) = / drT(r) = 1 + e 2 . (B6) 
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2. Velocity autocorrelation function of very hard particles 

Similarly to the Maxwell gas case, in order to determine the autocorrelation 
function T(r), we first have to solve kinetic equation fl32l) satisfied by the velocity 
distribution with the initial condition G(w; 0) = G in (w) = (w— < w > st )G st (w), and 
afterwards we have to evaluate the first moment of that solution 



J dw w G(w; r) = T(r) . 



Since norm Nq(t) vanishes, kinetic equation fl32l) becomes 

— + e— ) G(w; r) = [N 2 (r) - 2wN 1 (t)]$(w) - (w 2 + l)G(w; r) . (B7) 



^dr dw ^ 

That equation can be rewritten in Laplace world as 



eG(w; z) = / du exp {[S(u; z) — S(w; z)]/e} 

J — oo 

x{(u-<w> st )G st (u) + [N 2 (z)-2N 1 (z)uMu)} . (B8) 

The zeroth and the first moments of equation ( 1B8I) provide the system of equations 

= Joi + iV 2 A 00 - 2iMoi (B9) 
eNi = Jn + N 2 A W - 2NxA n 
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where 



dw / du exp {[S(u; z) — S(w; z)]/e} (u— < w > st )G st (u) (BIO) 

J — oo 

and 

Jn = / du» / du exp {[S(u; z) — S(w; z)]/e} w (u— < w > st )G st (u) . (Bll) 



The first moment N\(z) = T(z) is found to be 

f r \ = A 00 (z)Ju(z) - A 10 (z)Jqi(z) 

{ ] eA 00 (z) + 2[A 00 (z)A n (z) - A 10 (z)A 01 (z)} ' 1 ] 

where the shorthand notation Aji(z) = A,/(0;2:) has been used. The value of T(z) 
at z = yields the diffusion coefficient -D(e). We find here the same formula as that 
derived from the analysis of the hydrodynamic pole (see expression (1751)). 



Appendix C: Useful integral expressions for functions arising in the case of 
very hard particles 

In the integral representation (1A6j) of function Aji(k; z), it is useful to make the 
variable change u = w — ey. This leads to a double integral of the form 

oo 



dw / dylji(w,y) . 
Jo 

Thanks to the fast decay of integrand Iji(w, y) in any direction of plane (w, y), the 
integrals upon w and y can be exchanged. Then in the integral upon w, we make 
the variable change w — > rj with 

v = (i + 2y) 1/2 (w- y{1 + y) + ; ky , 

' K } V 1 + 2?/ e(l + 2j/) 

This provides 

f /' . 1 



Mk]z) = 7^l dy \ 



x exp I — (z + l)y — e 2 



7/ 3 (2 + t/) feV zfce?/ 2 (l + ?/) 
6(1 + 2y) 2(1 + 2?/) l + 2j/ 

t rf ( 7] ey(l+y) iky x y 

x / dr/exp(-— ) - — — - + 



2 , \l + 2y) 1 l 2 l + 2y l + 2y 

x( (CI) 

Vl + 2y) 1 /2 i + 2y l + 2yy v ; 
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The integral upon r\ can be easily performed, thanks to the simple dependence of the 
corresponding integrand with respect to 77, namely a Gaussian times a polynomial. 
The result is a combination of algebraic functions of y with coefficients which reduce 
to positive integer powers of k. Thus, the remaining integral upon y does converge 
for any complex value of k and z, thanks to the presence of factor 



exp 1 — e 2 



y 3 (2 + y) 



6(1 + 2y), 

which ensures a fast integrable decay of the integrand when y — > 00. That fast decay 
guarantees that Aji(k; z) is an entire function of both complex variables k and z. 

Integral representation ( 1C1I) can be specified to k = 0, j = 0, 1 and I = 0, 1. This 
provides useful expressions for functions A 00 (z), A w (z), A m (z) and A\\(z) which 
are analogous to formula (IMil for A(z). That formula is derived as follows. First, 
we compute dA/dz from expression flAlOl) specified to k = in terms of functions 
Aji(z) and of their derivatives with respecto z. Using the integral representations 
(IA6j) . such derivatives are then expressed in terms of the Aji(z) 7 s by combining 
differentiation under the integral sign and integration by parts. This allows us to 
infer that A(z) is the solution of the first order differential equation 

^_A = 3^-^„„, (C2) 
dz dz 

with the boundary condition at infinity A(z) — > when z = x — > 00. A straightfor- 
ward application of the constant-variation method leads to 



A(z) = 3 eA 00 (z 



dz' exp(-z')A 00 (z') . (C3) 



Eventually, we use the above integral representation of Aqo(z') into expression (1C3f) . 
and we exchange integrals upon z' and y thanks to absolute convergence. Since the 
dependence in z' reduces to simple exponential factor exp(— z'{\ + y)), the integral 
upon z' is readily done, and this eventually leads to formula ( |54|) . 
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